/*
    <one line to give the program's name and a brief idea of what it does.>
    Copyright (C) <year>  <name of author>

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.

*/

#include "uqBayesianScalarModel.h"

namespace atomism
{
uqBayesianScalarModel* ThisInvPb;
using namespace std;


//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------
	    
double likelihoodRoutine_(const uqGslVectorClass& paramValues,const uqGslVectorClass* paramDirection, const void*  Data,
			 uqGslVectorClass*  gradVector,uqGslMatrixClass*  hessianMatrix,uqGslVectorClass* hessianEffect)
{

 double result=ThisInvPb->LikelihoodRoutine(paramValues,paramDirection,Data,gradVector,hessianMatrix,hessianEffect);
  
 return(result);
 }
  


bool uqBayesianScalarModel::isBayesianScalarModelRegisteredInPython=0;
msRegistrar uqBayesianScalarModel::Registrar( "uqBayesianScalarModel" , uqBayesianScalarModel::createInstance );

//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

std::ostream& operator<<(std::ostream& out,const uqBayesianScalarModel& bayes) {
  
 return( bayes.print(out) );
 }

//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

void uqBayesianScalarModel::registryInPython(){
 
 #if USE_PYTHON
 msPhysicalInterface::registryInPython();

 if(  ! uqBayesianScalarModel::isBayesianScalarModelRegisteredInPython ) {
   
    using namespace boost::python;
	  
    class_<uqBayesianScalarModel,bases<msPhysicalInterface>,boost::shared_ptr<uqBayesianScalarModel> >
          ("uqBayesianScalarModel",  "Bayesian inferance of scalar model",no_init	  
	    ) 
    .def( "New",&uqBayesianScalarModel::New ,"Return a shared pointer to a new object" )
	.staticmethod("New")
    .def( "setDatabase",&uqBayesianScalarModel::setDatabase ,"Set the database. arg2: database" )
    .def( "setParameter",&uqBayesianScalarModel::setParameter ,
	  "Set or add a parameter, this parameter has to exist as constant in the model. arg2: name, arg3: min, arg4: max" )
    .def( "setModel",&uqBayesianScalarModel::setModel ,"Set the model. arg2: model" ) 
    .def( "setSampleFromJointDistribution",&uqBayesianScalarModel::setSampleFromJointDistribution ,
	  "Set the parameters from a sample of the joint distribution. arg2: indice of the sample" )
    .def( "infer",&uqBayesianScalarModel::infer ,"Infer the model. arg2: if true ask before computation" )
    .def( "getUqParamaters",&uqBayesianScalarModel::getUqParamaters ,"Return the parameter map" )
    .def( "getPosteriorPDF",&uqBayesianScalarModel::getPosteriorPDF ,"Return the posterior PDF" );
    uqBayesianScalarModel::isBayesianScalarModelRegisteredInPython = 1 ;
    }
#endif
}

//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

uqBayesianScalarModel::uqBayesianScalarModel() {
 
 int initialized; 
 MPI_Initialized(&initialized);
 if (!initialized) MPI_Init(NULL, NULL);
 constructVar("uqBayesianScalarModel","BayesianScalarModel","bayesian for scalar model");
 }
    
//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

boost::shared_ptr<msTreeMapper> uqBayesianScalarModel::setModel(boost::shared_ptr<msScalarFunction> model){
 
  if(!model->hasParent()) model->setAffiliation(mySharedPtr());
  msTreeMapper::update(Model,model);
  PosteriorPDF->setModel(Model.getSharedPtr());
  return mySharedPtr();
}

//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

boost::shared_ptr<msTreeMapper> uqBayesianScalarModel::setDatabase(boost::shared_ptr<msDiracs> database){
  
  if(!database->hasParent()) database->setAffiliation(mySharedPtr());

  msTreeMapper::update(Database,database);
  return mySharedPtr();
}
//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

boost::shared_ptr<msTreeMapper> uqBayesianScalarModel::setParameter(string name,double min,double max) {
  
 if(Model==boost::shared_ptr<msScalarFunction>() )
   BOOST_THROW_EXCEPTION(msError("you need to define the model before the parameters",
				 "void uqBayesianScalarModel::setParameter(size_t i,double min,double max)",
				 getFullId()));

 try{ double& ref = Model->getConstant(name);
      if( (ref>max) || (ref<min) )
        BOOST_THROW_EXCEPTION(msError("The nominal value has to be between min and max values",
				      "void uqBayesianScalarModel::setParameter(size_t i,double min,double max)",
				       getFullId()));
      boost::shared_ptr<msScalarVariable> param=msScalarVariable::New("");
      param->set(ref,min,max,0,0);
      param->Default=ref;
      param->setId(name);
      UqParameters->addVar(param);
  }
  catch(msError& e){
    e.addContext("the constant "+name+" does not exist (boost::shared_ptr<msTreeMapper> uqBayesianScalarModel::setParameter(size_t i,double min,double max)");    
  }
 return mySharedPtr();
 }
  
//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

boost::shared_ptr<msTreeMapper> uqBayesianScalarModel::setSampleFromJointDistribution(size_t j){
  
 try{
   PosteriorPDF->setSample(UqParameters.getSharedPtr(),j); 
 }
 catch(msError& e) { 
   e.addContext("can not set the sample from joint distribution (boost::shared_ptr<msTreeMapper> uqBayesianScalarModel::setSampleFromJointDistribution(size_t j))");
 }
 updateModelParameters();
 return mySharedPtr();
 }
   
//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

boost::shared_ptr<msTreeMapper> uqBayesianScalarModel::updateModelParameters(){
  
 for( size_t i=0;i<UqParameters->noOfDim();++i)
     Model->setConstant( (*UqParameters)[i]->getId(), (*UqParameters)[i]->Value);
 return mySharedPtr();
 }
 
//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

double uqBayesianScalarModel::LikelihoodRoutine(const uqGslVectorClass& paramValues,
						const uqGslVectorClass* paramDirection,
						const void* functionDataPtr,
						uqGslVectorClass* gradVector,
						uqGslMatrixClass* hessianMatrix, 
						uqGslVectorClass* hessianEffect) {
  
  if(paramValues.sizeGlobal()!=UqParameters->noOfDim() + 1)
     BOOST_THROW_EXCEPTION(msError("The 'paramValues' vector's size is not coherent with the number of parameters.",
				   "double uqBayesianScalarModel::LikelihoodRoutine(...)",
				   getFullId()));  

  double sigma2 = paramValues[UqParameters->noOfDim()];   
  
  for( size_t i=0;i<UqParameters->noOfDim();++i)
     (*UqParameters)[i]->Value=paramValues[i]* (*UqParameters)[i]->Default;
  
  updateModelParameters();

  boost::shared_ptr<msSamplerExplicit> DbSampler=Database->getSampler();

  double likeli = 0;
  DbSampler->begin();
  while(DbSampler->getNextPoint(Model->getCoordinates())){
    
     double error  = fabs(Database->evaluate()-Model->evaluate());
     likeli = likeli - 0.5*log( 2.0 * csts::pi * sigma2) - error*error/sigma2/2.0;
  }
  
  if( ( CurrentIteration % DisplayStep == 0 ) || ( likeli>LikelihoodMax ) ){

    stringstream out;    
    if(likeli>LikelihoodMax){ out<<"New maximum of likelihood"<<endl;LikelihoodMax=likeli; }

    for(size_t i=0;i<UqParameters->noOfDim();i++)  out<<setw(15)<<(*UqParameters)[i]->Value;
    out<<setw(15)<<sigma2<<setw(15)<<likeli;
    LOGGER_WRITE( msLogger::INFO, out.str() );
  }
 
  CurrentIteration++;
  return likeli;
 };

//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

void uqBayesianScalarModel::infer() {

  LOGGER_ENTER_FUNCTION("void uqBayesianScalarModel::infer(bool askok)",getFullId())

  if( (Database != boost::shared_ptr<msDiracs>()) &&
      (Database->getCoordinates() != Model->getCoordinates() ) )
     BOOST_THROW_EXCEPTION(msError("The model and the database have to share the same coordinates",
				   "void uqBayesianScalarModel::infer()",
				   getFullId())); 
  if( UqParameters->noOfDim() == 0 )
     BOOST_THROW_EXCEPTION(msError("There is no parameters defined, please define them by calling 'setParameter'.",
				   "void uqBayesianScalarModel::infer()",
				   getFullId()));

  boost::shared_ptr<msSamplerExplicit> DbSampler=Database->getSampler();
  DbSampler->begin();
  if(! DbSampler->getNextPoint(Model->getCoordinates()) )  
    BOOST_THROW_EXCEPTION(msError("There is no data loaded in the database, please call 'loadDatabase'",
				   "void uqBayesianScalarModel::testLikelihood()",
				   getFullId())); 
  writeInputQueso();
  
  Environment = boost::shared_ptr<uqFullEnvironmentClass>(
                new uqFullEnvironmentClass(MPI_COMM_WORLD,"QuesoInput.inp","",NULL)
                 );
  CurrentIteration=0;
  uqVectorSpaceClass<uqGslVectorClass,uqGslMatrixClass>  
		paramSpace(*Environment, "param_", UqParameters->noOfDim() + 1, NULL);
  uqGslVectorClass paramMins(paramSpace.zeroVector());
  uqGslVectorClass paramMaxs(paramSpace.zeroVector());

  for(size_t i=0;i<UqParameters->noOfDim();i++ ){
    
    paramMins[i]=(*UqParameters)[i]->Min/(*UqParameters)[i]->Default;
    paramMaxs[i]=(*UqParameters)[i]->Max/(*UqParameters)[i]->Default;
  }
  paramMins[UqParameters->noOfDim()] = 0;
  paramMaxs[UqParameters->noOfDim()] = 1;

  uqBoxSubsetClass<uqGslVectorClass,uqGslMatrixClass>
  paramDomain("param_",paramSpace,paramMins,paramMaxs);
  uqGslVectorClass meanVector(paramSpace.zeroVector());

  uqGenericScalarFunctionClass<uqGslVectorClass,uqGslMatrixClass>
    			likelihoodFunctionObj(	"like_",
                          			paramDomain,
                          			likelihoodRoutine_,
                          			(void*) this,
                          			true); 
  
  uqUniformVectorRVClass<uqGslVectorClass,uqGslMatrixClass> 
			priorRv("prior_", paramDomain);
  uqGenericVectorRVClass<uqGslVectorClass,uqGslMatrixClass> 
			postRv("post_", paramSpace);
  uqStatisticalInverseProblemClass<uqGslVectorClass,uqGslMatrixClass>
   		  	ip("", NULL, priorRv, likelihoodFunctionObj, postRv);
   
  stringstream out; int l=15; 
  for(int i=0;i<UqParameters->noOfDim();i++) {
     
    if( (*UqParameters)[i]->getId().length()>10) 
         out<<setw(15)<<(*UqParameters)[i]->getId();
  }
  LOGGER_WRITE( msLogger::INFO, out.str() );

  ThisInvPb=this;
  LikelihoodMax = -100000000000;
  int lChain; double ratio;
  getParameters()->queryDoubleValue("DisplayStepRatio", ratio, *getUnits());
  getParameters()->queryIntValue(   "ChainLength"     , lChain);
  DisplayStep = lChain * ratio;

  ip.solveWithBayesMLSampling();

 
  uqGslVectorClass tmpVec(postRv.imageSet().vectorSpace().zeroVector());
  PosteriorPDF->clear();
  for (unsigned int i = 0; i < postRv.realizer().subPeriod(); ++i) {
    postRv.realizer().realization(tmpVec);
    vector<double>LastChain_i;
    for(int j = 0; j<UqParameters->noOfDim() ; j++)
      LastChain_i.push_back( tmpVec[j]*(*UqParameters)[j]->Default );
     
    PosteriorPDF->addSample(LastChain_i);
  }
    
  
 LOGGER_EXIT_FUNCTION();
 }

//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

std::ostream& uqBayesianScalarModel::print(std::ostream& out) const
{
 msPhysicalInterface::print(out);
 output::printHeader(out,"Bayesian inferance");

 out<<"Number of parameters: "<<UqParameters->noOfDim()<<endl;
 out<<setw(15)<<"name"<<setw(15)<<"min"<<setw(15)<<"max"<<setw(15)<<"nominal"<<setw(25)<<"value in model"<<endl;

 for(size_t i=0;i<UqParameters->noOfDim();++i)
    out<<setw(15)<<(*UqParameters)[i]->getId()<<setw(15)<<(*UqParameters)[i]->Min<<setw(15)<<(*UqParameters)[i]->Max
       <<setw(15)<<(*UqParameters)[i]->Default<<setw(25)<<Model->getConstant((*UqParameters)[i]->getId())<<endl;

 out<<endl<<"The database contains "<<Database->noOfValues()<<" values"<<setprecision(8)<<endl;

 msSamplerExplicit& sampler = *(Database->getSampler());
 msGeneralizedCoordinates& gcoor = *(Database->getCoordinates());

 sampler.begin();
 for(size_t i=0;i<gcoor.noOfDim();++i)  out<<setw(15)<<gcoor[i]->getId();
 out<<setw(15)<<"Value"<<setw(15)<<"Model"<<endl; 

 while(sampler.getNextPoint( Database->getCoordinates() )){

    if( (sampler.noOfScanned()<10) || (sampler.noOfScanned() > Database->noOfValues()-10) ){

       for(size_t i=0;i<gcoor.noOfDim();++i)  out<<setw(15)<<gcoor[i]->Value;
       out<<setw(15)<<Database->evaluate()<<setw(15)<<Model->evaluate()<<endl;
       }
    if( (Database->noOfValues()>20) &&(sampler.noOfScanned() == 10) ) out<<setw(15)<<".\n"<<setw(15)<<".\n"<<setw(15)<<".\n"<<endl;
    }        
 
 return out;
 }
 
//-------------------------------------------------------------------------------------------------
//-------------------------------------------------------------------------------------------------

void uqBayesianScalarModel::writeInputQueso() {
  
 int L,l;
 getParameters()->queryIntValue("ChainLength",l); 
 getParameters()->queryIntValue("ChainLength",L); 

 std::ofstream input("QuesoInput.inp");
 input<<"###############################################\n"
 <<"## UQ Environment\n"
 <<"################################################\n"
 <<"#env_help                = anything\n"
 <<"env_numSubEnvironments   = 1\n"
 <<"env_subDisplayFileName   = display\n"
 <<"env_subDisplayAllowAll   = 0\n"
 <<"env_subDisplayAllowedSet = 0\n"
 <<"env_displayVerbosity     = 0\n"
 <<"env_syncVerbosity        = 0\n"
 <<"env_seed                 = 0\n"
 <<"################################################\n"
 <<"## Statistical inverse problem (ip)\n"
 <<"################################################\n"
 <<"ip_help                 = anything\n"
 <<"ip_computeSolution      = 1\n"
 <<"ip_dataOutputFileName   = sipOutput\n"
 <<"ip_dataOutputAllowedSet = 0\n"
 <<"################################################\n"
 <<"## 'ip_': information for Multilevel algorithm\n"
 <<"################################################\n"
 <<"ip_ml_help                 = anything\n"
 <<"ip_ml_dataOutputFileName   = sipOutput_ml\n"
 <<"ip_ml_dataOutputAllowedSet = 0 1\n"
 <<"################################################\n"
 <<"## All levels, unless otherwise specified\n"
 <<"################################################\n"
 <<"ip_ml_default_help=xyz\n"
 <<"ip_ml_default_rawChain_size ="<<l<<"\n"
 <<"ip_ml_default_rawChain_dataOutputFileName = rawChain_ml\n"
 <<"ip_ml_0_rawChain_dataOutputFileName = rawChain_ml\n"
 <<"ip_ml_default_putOutOfBoundsInChain = 1\n"
 <<"ip_ml_default_scaleCovMatrix                      = 1\n"
 <<"ip_ml_default_dr_maxNumExtraStages                = 0\n"
 <<"################################################\n"
 <<"## Last level (level 4 for 'example_gsl' executable)\n"
 <<"################################################\n"
 <<"ip_ml_last_dataOutputFileName                  = sipOutput_ml\n"
 <<"ip_ml_last_dataOutputAllowedSet                = 0  1\n"
 <<"ip_ml_last_rawChain_size                       = "<<L<<"\n"
 <<"ip_ml_last_rawChain_computeStats               = 1\n"
 <<"ip_ml_last_rawChain_dataOutputFileName         = rawChain_ml\n"
 <<"ip_ml_last_rawChain_stats_kde_compute          = 1\n"
 //<<"ip_ml_last_rawChain_stats_kde_numEvalPositions = 1000\n"
 <<"ip_ml_last_rawChain_stats_covMatrix_compute    = 1\n"
 <<"ip_ml_last_rawChain_stats_corrMatrix_compute   = 1\n"
 <<"ip_ml_last_putOutOfBoundsInChain = 1\n"
 <<"ip_ml_last_scaleCovMatrix                      = 1\n"
 <<"ip_ml_last_dr_maxNumExtraStages                = 0\n";
 }

}


/*
  int     D_number = 100; // The number of the data points
  
  std::ofstream outdata;

  uqVectorSpaceClass<uqGslVectorClass,uqGslMatrixClass> CovSpace(env, "", D_number, NULL);
  uqGslMatrixClass* mat1=CovSpace.newMatrix();

  double solution_vector[D_number];
  double ExpData_vector[D_number];

  // definition of the covariance matrix of mat1
  for (int j = 0; j < D_number; ++j) {
     for (int i = 0; i < D_number; ++i) {
        if(i==j){
           if(i <  D_number) (*mat1)(i,j) =//sigma^2  cal_sigmaML_Sq; 
        } 
        if(i > j) {
           if(Correlation_UQ_model==1) { 
              if(i <   D_number  and j <  D_number)(*mat1)(i,j) = cal_sigmaML_Sq * exp(-pow(fabs(Exp_data_mol_frac_N1[i]-Exp_data_mol_frac_N1[j])/cal_ML_Clength,cal_ML_power)); 
              if(i >=  D_number  and j >= D_number)(*mat1)(i,j) = cal_sigmaXN_Sq * exp(-pow(fabs(Exp_data_mol_frac_N1[i-D_number]-Exp_data_mol_frac_N1[j-D_number])/cal_XN_Clength,cal_XN_power));
           }
           if(Correlation_UQ_model==2) {
              if(i - D_number  ==  j) (*mat1)(i,j) = cal_sigmaML * cal_sigmaXN * cal_MLXN_Correlation; 
           }
           if(Correlation_UQ_model==3) {
              if(i <   D_number and j <  D_number)(*mat1)(i,j) = cal_sigmaML_Sq * exp(-pow(fabs(Exp_data_mol_frac_N1[i]-Exp_data_mol_frac_N1[j])/cal_ML_Clength,cal_ML_power));              
              if(i >=  D_number and j >= D_number)(*mat1)(i,j) = cal_sigmaXN_Sq * exp(-pow(fabs(Exp_data_mol_frac_N1[i-D_number]-Exp_data_mol_frac_N1[j-D_number])/cal_XN_Clength,cal_XN_power));
              if(i >=  D_number and j <  D_number)(*mat1)(i,j) = cal_sigmaML * cal_sigmaXN * cal_MLXN_Correlation * exp(-pow(fabs(Exp_data_mol_frac_N1[i-D_number]-Exp_data_mol_frac_N1[j])/cal_MLXN_Clength,cal_MLXN_power));  
           }
           (*mat1)(j,i) = (*mat1)(i,j);
        }
     }
  }

  // obtain the inverse matrix of mat1 
  uqGslMatrixClass matI(*mat1);
  matI = mat1->inverse();



  for(unsigned int i = 0; i < D_number; ++i) {
        solution_vector[i] = error [i]// your model output - data; 
        ExpData_vector[i]  =// data;
  }   // i = Exp, index

  resultValue = 0.0;
  double tmp;
  double tmp_vec[D_number];
  double determinant = 0.0;
  determinant = mat1->determinant();
  for (unsigned int i = 0; i < D_number; ++i) {
     tmp_vec[i]=0.0;
     tmp = 0;
     for (unsigned int j = 0; j < D_number; ++j) {
        tmp_vec[i] += 0.5 * matI(j,i) * solution_vector[j];
     }
     if(cal_UQ_model==0) resultValue += -0.5 * log(2.0*pi)- tmp_vec[i]*solution_vector[i];
  }

// Compute log_det_K = log(det(mat1)) (use Cholesky to make >> this easy)
  int ierr = mat1->chol();
  if( ierr != 0 ) printf("Cholesky decomposition failed!\n"); //after cholesky transformation: det(mat1)=2.*sum_i mat[i][i]
  double log_det_K = 0.0;
  for( int ii=0; ii<2 * D_number; ii++ ) log_det_K += 2.0*std::log((*mat1)(ii,ii));

  resultValue   += -0.5 * log_det_K;

  if(ierr != 0) {resultValue =-100000;}

  delete mat1;
  

  return resultValue;
}
*/
